vcs <- rio::import("data_complete_all_anonym_allZ.xlsx")
var_label(vcs$dominance) <- "Dominance"
var_label(vcs$neuro) <- "Neuroticism"
var_label(vcs$agree) <- "Agreeableness"
var_label(vcs$extra) <- "Extraversion"
var_label(vcs$openn) <- "Openness"
var_label(vcs$consc) <- "Conscientiousness"
var_label(vcs$soi_full) <- "Unrestricted sociosexuality"
var_label(vcs$f0) <- "Voice pitch"
var_label(vcs$pf) <- "Formants"
vcs$sex_c <- vcs$sex
vcs$sex <- factor(if_else(vcs$sex == 1, "male", "female"))
contrasts(vcs$sex) <- contr.helmert(2)
var_label(vcs$sex) <- "Sex"
set.seed(1)
var_label(vcs$age) <- "Age"
vcs <- vcs %>%
mutate(
age = if_else(dataset == 9, 20, age)/10,
age_se = if_else(dataset == 9, 3, 0.5)/10
)
warmup <- 2000
iter <- warmup + 2000
chains <- 4
control <- list(adapt_delta = 0.99)
priors <- c(
prior(normal(0, 3), class = b)
)
variable_labels <- c("Intercept"= "Intercept", "sex1" = "Sex [male]", "bsp_meageage_se" = "Age±SE", "f0" = "Voice pitch (f0)", "pf" = "Formants (f1-f4)", "age" = "Age")
effect_labels <- c("b_Intercept"= "Intercept", "b_sex1" = "Sex [male]", "bsp_meageage_se" = "Age±SE", "b_f0" = "Voice pitch (f0)", "b_pf" = "Formants (f1-f4)", "b_age" = "Age")
codebook::md_pattern(vcs %>% select(dominance, sex, age, f0, pf, dataset))
## # A tibble: 4 x 7
## description age pf f0 dominance var_miss n_miss
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Missing values per variable 4 4 7 1256 1271 1271
## 2 Missing values in 1 variables 1 1 1 0 1 1245
## 3 Missing values in 0 variables 1 1 1 1 0 985
## 4 3 other, less frequent patterns 2 2 1 0 7 11
vcs <- vcs %>% filter(!is.na(f0), !is.na(pf), !is.na(age))
Participants with lower voice pitch will have a more dominant personality.
Interpretation: The data are are consistent with substantial linear negative relationships between f0 and dominance. The 89% HDI for f0 falls entirely outside the ROPE, but the HDI for pf falls almost entirely within it. This evidence is consistent with a non-negligible association, where people with deeper voices are more dominant (after adjusting for age and gender). Further research is needed to verify if the association with pf is truly negligible in size.
ggplot(vcs, aes(f0, dominance)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing missing values (geom_point).
ggplot(vcs, aes(pf, dominance)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing missing values (geom_point).
ggplot(vcs, aes(age*10, dominance)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
scale_x_continuous("Age") +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing non-finite values (stat_smooth).
## Warning: Removed 1245 rows containing missing values (geom_point).
Decisions upon visual diagnosis: Visual diagnostic show no clear sign of nonlinearities. We will fit linear effects for f0, pf, and age. It is not necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is not included in these analyses. Visually, there could be an interaction between sex and f0.
h1_simple <- brm(dominance ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
iter = iter, warmup = warmup, chains = chains, cores = chains,
prior = priors,
control = control, save_mevars = TRUE,
file = "models/dominance/h1_simple") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/dominance/h1_simple.rds'
## Automatically saving the model object in 'models/dominance/h1_simple.rds'
## Automatically saving the model object in 'models/dominance/h1_simple.rds'
h1_after_diagnosis <- brm(dominance ~ f0*sex + pf + sex + age + (1 | dataset), data = vcs,
iter = iter, warmup = warmup, cores = chains, chains = chains,
prior = priors,
control = control, save_mevars = TRUE,
file = "models/dominance/h1_nonlinear") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/dominance/h1_nonlinear.rds'
## Automatically saving the model object in 'models/dominance/h1_nonlinear.rds'
## Automatically saving the model object in 'models/dominance/h1_nonlinear.rds'
# conditional_smooths(h1_after_diagnosis)
compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
## elpd_diff se_diff
## h1_simple 0.0 0.0
## h1_after_diagnosis -1.0 0.2
h1_sex_mod_dataset_varying <- brm(dominance ~ f0 + sex + age +
(1 + f0 || dataset), data = vcs,
prior = priors,
iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
control = control, save_mevars = TRUE,
file = "models/dominance/h1_sex_mod_dataset_varying") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/dominance/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/dominance/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/dominance/h1_sex_mod_dataset_varying.rds'
h1_after_diagnosis_melsm <- brm(
bf(dominance ~ f0 + sex + age + (1 | dataset),
sigma ~ f0 + sex + age + (1 | dataset)), data = vcs,
prior = priors,
iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
control = control, save_mevars = TRUE,
file = "models/dominance/h1_after_diagnosis_melsm") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/dominance/h1_after_diagnosis_melsm.rds'
## Automatically saving the model object in 'models/dominance/h1_after_diagnosis_melsm.rds'
## Automatically saving the model object in 'models/dominance/h1_after_diagnosis_melsm.rds'
compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying,
h1_after_diagnosis_melsm)
compare_models
## elpd_diff se_diff
## h1_after_diagnosis_melsm 0.0 0.0
## h1_simple -0.2 4.0
## h1_after_diagnosis -1.2 4.0
## h1_sex_mod_dataset_varying -2.9 4.4
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying,
h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>%
map(~ { x <- as_tibble(bayes_R2(.))
x$loo_R2 <- loo_R2(.)
x
}) %>%
bind_rows(.id="model") %>%
arrange(desc(loo_R2))
sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
pred.labels = variable_labels)
| Â | Dominance | |
|---|---|---|
| Predictors | Estimates | CI (95%) |
| Intercept | -0.08 | -0.80;0.60 |
| Voice pitch (f0) | -0.27 | -0.46;-0.08 |
| Formants (f1-f4) | -0.02 | -0.20;0.16 |
| Sex [male] | -0.31 | -0.52;-0.08 |
| Age±SE | 0.00 | -0.01;0.02 |
| Random Effects | ||
| σ2 | 0.05 | |
| τ00 | 0.96 | |
| ICC | 0.05 | |
| N dataset | 4 | |
| Observations | 985 | |
| Marginal R2 / Conditional R2 | 0.069 / 0.069 | |
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
## Possible multicollinearity between b_sex1 and b_f0 (r = 0.76). This might lead to inappropriate results. See 'Details' in '?equivalence_test'.
equiv_test_table <- equiv_test %>%
mutate(Parameter = recode(Parameter, !!!effect_labels)) %>%
mutate(
ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>%
select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
| Parameter | H0 | inside ROPE | 89% HDI |
|---|---|---|---|
| Intercept | Undecided | 32% | [-0.60; 0.41] |
| Voice pitch (f0) | Rejected | 0% | [-0.42;-0.11] |
| Formants (f1-f4) | Undecided | 81% | [-0.16; 0.13] |
| Sex [male] | Rejected | 0% | [-0.50;-0.14] |
| Age±SE | Accepted | 100% | [-0.01; 0.01] |
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Picking joint bandwidth of 0.00905
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).
Model details
summary(h1_simple)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 5 divergent transitions after warmup. Increasing adapt_delta
## above 0.99 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: dominance ~ f0 + pf + sex + me(age, age_se) + (1 | dataset)
## Data: vcs (Number of observations: 985)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~dataset (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.52 0.43 0.16 1.66 1.00 1854 2414
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.08 0.37 -0.80 0.60 1.00 1470 1979
## f0 -0.27 0.10 -0.46 -0.08 1.00 4886 5623
## pf -0.02 0.09 -0.20 0.16 1.00 7512 6543
## sex1 -0.31 0.11 -0.52 -0.08 1.00 4367 5610
## meageage_se 0.00 0.01 -0.01 0.02 1.00 7099 5994
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.97 0.02 0.93 1.01 1.00 8321 5117
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta
## above 0.99 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: dominance ~ f0 * sex + pf + sex + me(age, age_se) + (1 | dataset)
## Data: vcs (Number of observations: 985)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~dataset (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.53 0.43 0.16 1.63 1.00 2619 3493
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.05 0.37 -0.79 0.67 1.00 3111 3051
## f0 -0.26 0.10 -0.45 -0.07 1.00 6625 5759
## sex1 -0.30 0.12 -0.52 -0.07 1.00 5716 5773
## pf -0.02 0.09 -0.19 0.16 1.00 9331 6207
## f0:sex1 0.02 0.10 -0.17 0.21 1.00 8589 6232
## meageage_se 0.00 0.01 -0.01 0.02 1.00 9747 6170
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.97 0.02 0.93 1.02 1.00 10519 5550
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_after_diagnosis_melsm)
# conditional_smooths(h1_after_diagnosis_melsm)
summary(h1_sex_mod_dataset_varying)
## Warning: There were 4 divergent transitions after warmup. Increasing adapt_delta
## above 0.99 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: dominance ~ f0 + pf + f0 * sex + me(age, age_se) + (1 + f0 * sex || dataset)
## Data: vcs (Number of observations: 985)
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~dataset (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.50 0.47 0.03 1.78 1.00 3660 3397
## sd(f0) 0.26 0.35 0.01 1.20 1.00 3195 3060
## sd(sex1) 0.34 0.43 0.01 1.51 1.00 2747 3104
## sd(f0:sex1) 0.37 0.39 0.02 1.41 1.00 2894 3588
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.09 0.39 -0.84 0.68 1.00 5857 4642
## f0 -0.25 0.26 -0.71 0.20 1.00 3824 2480
## pf -0.03 0.09 -0.21 0.16 1.00 11668 8835
## sex1 -0.29 0.30 -0.95 0.29 1.00 3959 2492
## f0:sex1 -0.00 0.28 -0.56 0.54 1.00 4020 3426
## meageage_se 0.00 0.01 -0.01 0.02 1.00 11857 9601
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.97 0.02 0.93 1.01 1.00 13009 8654
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis_melsm)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta
## above 0.99 may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup
## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: dominance ~ f0 * sex + pf + sex + me(age, age_se) + (1 | dataset)
## sigma ~ f0 * sex + pf + sex + me(age, age_se) + (1 | dataset)
## Data: vcs (Number of observations: 985)
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~dataset (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.52 0.39 0.17 1.58 1.00 3584 4787
## sd(sigma_Intercept) 0.15 0.17 0.01 0.62 1.00 2706 3542
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.07 0.35 -0.79 0.60 1.00 3595 3307
## sigma_Intercept 0.04 0.19 -0.34 0.38 1.00 4339 5652
## f0 -0.29 0.10 -0.48 -0.09 1.00 7558 8924
## sex1 -0.30 0.12 -0.54 -0.07 1.00 6195 7846
## pf 0.00 0.09 -0.16 0.18 1.00 11479 9318
## f0:sex1 0.04 0.10 -0.15 0.23 1.00 12251 9131
## sigma_f0 0.10 0.07 -0.04 0.25 1.00 7277 8700
## sigma_sex1 0.13 0.09 -0.04 0.30 1.00 6632 8402
## sigma_pf 0.03 0.07 -0.11 0.16 1.00 11282 9463
## sigma_f0:sex1 -0.01 0.07 -0.15 0.13 1.00 10718 8125
## meageage_se 0.00 0.01 -0.01 0.02 1.00 11342 9225
## sigma_meageage_se -0.00 0.01 -0.01 0.01 1.00 4908 7372
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
| model | Estimate | Est.Error | Q2.5 | Q97.5 | loo_R2 |
|---|---|---|---|---|---|
| h1_simple | 0.0699913 | 0.0147882 | 0.0427074 | 0.1001835 | 0.0529997 |
| h1_after_diagnosis | 0.0705955 | 0.0149998 | 0.0434176 | 0.1014367 | 0.0511312 |
| h1_after_diagnosis_melsm | 0.0723500 | 0.0145552 | 0.0455458 | 0.1020524 | 0.0509098 |
| h1_sex_mod_dataset_varying | 0.0754052 | 0.0149069 | 0.0474739 | 0.1061981 | 0.0476857 |
Participants with lower voice pitch will score higher on agreeableness.
Interpretation: The credible intervals for f0 and pf include zero, but the 89% HDI does not fall entirely within the ROPE. This evidence is consistent with a negligible association between f0/pf and agreeableness. Further research is needed to verify if the associations are truly negligible in size.
ggplot(vcs, aes(f0, agree)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).
ggplot(vcs, aes(pf, agree)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).
ggplot(vcs, aes(age*10, agree)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
scale_x_continuous("Age") +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).
Decisions upon visual diagnosis: Visual diagnostic shows signs of nonlinearity for f0 and pf, but no clear interactions by sex. We will fit splines for f0, and pf. It is necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is included in these analyses.
h1_simple <- brm(agree ~ f0 + pf + sex + me(age, age_se) + (1 | dataset), data = vcs,
iter = iter, warmup = warmup, chains = chains, cores = chains,
prior = priors,
control = control, save_mevars = TRUE,
file = "models/agree/h1_simple") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/agree/h1_simple.rds'
## Automatically saving the model object in 'models/agree/h1_simple.rds'
## Automatically saving the model object in 'models/agree/h1_simple.rds'
h1_after_diagnosis <- brm(agree ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset), data = vcs,
iter = iter, warmup = warmup, cores = chains, chains = chains,
prior = priors,
control = control, save_mevars = TRUE,
file = "models/agree/h1_nonlinear") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/agree/h1_nonlinear.rds'
## Automatically saving the model object in 'models/agree/h1_nonlinear.rds'
## Automatically saving the model object in 'models/agree/h1_nonlinear.rds'
conditional_smooths(h1_after_diagnosis)
compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
## elpd_diff se_diff
## h1_after_diagnosis 0.0 0.0
## h1_simple -3.8 3.1
h1_sex_mod_dataset_varying <- brm(agree ~ f0 + pf + f0 * sex + me(age, age_se) +
(1 + f0 + pf + sex || dataset), data = vcs,
prior = priors,
iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
control = control, save_mevars = TRUE,
file = "models/agree/h1_sex_mod_dataset_varying") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/agree/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/agree/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/agree/h1_sex_mod_dataset_varying.rds'
h1_after_diagnosis_melsm <- brm(
bf(agree ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset),
sigma ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset)), data = vcs,
prior = priors,
iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
control = control, save_mevars = TRUE,
file = "models/agree/h1_after_diagnosis_melsm") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/agree/h1_after_diagnosis_melsm.rds'
## Automatically saving the model object in 'models/agree/h1_after_diagnosis_melsm.rds'
## Automatically saving the model object in 'models/agree/h1_after_diagnosis_melsm.rds'
compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying,
h1_after_diagnosis_melsm)
compare_models
## elpd_diff se_diff
## h1_after_diagnosis_melsm 0.0 0.0
## h1_after_diagnosis -7.0 5.1
## h1_sex_mod_dataset_varying -7.5 6.2
## h1_simple -10.9 5.9
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying,
h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>%
map(~ { x <- as_tibble(bayes_R2(.))
x$loo_R2 <- loo_R2(.)
x
}) %>%
bind_rows(.id="model") %>%
arrange(desc(loo_R2))
sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
pred.labels = variable_labels)
| Â | Agreeableness | |
|---|---|---|
| Predictors | Estimates | CI (95%) |
| Intercept | 0.09 | -0.37;0.59 |
| Voice pitch (f0) | 0.01 | -0.13;0.16 |
| Formants (f1-f4) | 0.09 | -0.06;0.24 |
| Sex [male] | 0.02 | -0.15;0.20 |
| Age±SE | -0.01 | -0.02;0.00 |
| Random Effects | ||
| σ2 | 0.11 | |
| τ00 | 0.90 | |
| ICC | 0.11 | |
| N dataset | 7 | |
| Observations | 1434 | |
| Marginal R2 / Conditional R2 | 0.117 / 0.117 | |
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
## Possible multicollinearity between b_sex1 and b_f0 (r = 0.73). This might lead to inappropriate results. See 'Details' in '?equivalence_test'.
equiv_test_table <- equiv_test %>%
mutate(Parameter = recode(Parameter, !!!effect_labels)) %>%
mutate(
ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>%
select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
| Parameter | H0 | inside ROPE | 89% HDI |
|---|---|---|---|
| Intercept | Undecided | 35% | [-0.28; 0.47] |
| Voice pitch (f0) | Undecided | 91% | [-0.11; 0.13] |
| Formants (f1-f4) | Undecided | 58% | [-0.04; 0.20] |
| Sex [male] | Undecided | 80% | [-0.13; 0.16] |
| Age±SE | Accepted | 100% | [-0.01; 0.00] |
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Picking joint bandwidth of 0.00726
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).
Model details
summary(h1_simple)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: agree ~ f0 + pf + sex + me(age, age_se) + (1 | dataset)
## Data: vcs (Number of observations: 1434)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~dataset (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.48 0.21 0.25 0.99 1.00 1992 3446
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.09 0.25 -0.37 0.59 1.00 2017 3031
## f0 0.01 0.07 -0.13 0.16 1.00 6296 5690
## pf 0.09 0.08 -0.06 0.24 1.00 7233 5707
## sex1 0.02 0.09 -0.15 0.20 1.00 5467 5890
## meageage_se -0.01 0.01 -0.02 0.00 1.00 8114 6260
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.94 0.02 0.91 0.98 1.00 11815 6095
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: agree ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset)
## Data: vcs (Number of observations: 1434)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sf0_1) 1.05 0.65 0.23 2.73 1.00 3297 3480
## sds(spf_1) 0.65 0.50 0.03 1.91 1.00 3618 3082
##
## Group-Level Effects:
## ~dataset (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.48 0.20 0.24 1.01 1.00 2375 2765
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.11 0.24 -0.35 0.60 1.00 2712 3234
## sex1 0.13 0.13 -0.11 0.39 1.00 7526 6076
## sf0_1 -0.24 1.58 -3.54 2.84 1.00 6436 5755
## spf_1 0.60 1.26 -2.04 3.33 1.00 5442 4987
## meageage_se -0.01 0.01 -0.02 0.00 1.00 10169 6084
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.94 0.02 0.90 0.97 1.00 14100 6179
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_after_diagnosis_melsm)
# conditional_smooths(h1_after_diagnosis_melsm)
summary(h1_sex_mod_dataset_varying)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: agree ~ f0 + pf + f0 * sex + me(age, age_se) + (1 + f0 + pf + sex || dataset)
## Data: vcs (Number of observations: 1434)
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~dataset (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.53 0.23 0.26 1.14 1.00 3925 5371
## sd(f0) 0.15 0.12 0.01 0.44 1.00 3694 4711
## sd(pf) 0.14 0.12 0.01 0.44 1.00 4726 5476
## sd(sex1) 0.15 0.16 0.01 0.57 1.00 3137 4980
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.22 0.28 -0.34 0.77 1.00 3801 5048
## f0 0.04 0.11 -0.17 0.25 1.00 6903 5540
## pf 0.04 0.11 -0.18 0.24 1.00 7535 6008
## sex1 0.07 0.15 -0.23 0.35 1.00 5261 4665
## f0:sex1 0.13 0.08 -0.04 0.30 1.00 10675 8696
## meageage_se -0.01 0.01 -0.02 0.01 1.00 11214 8721
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.94 0.02 0.90 0.97 1.00 16491 8626
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis_melsm)
## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: agree ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset)
## sigma ~ s(f0) + s(pf) + sex + me(age, age_se) + (1 | dataset)
## Data: vcs (Number of observations: 1434)
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
## total post-warmup samples = 12000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sf0_1) 1.24 0.77 0.27 3.20 1.00 3949 5731
## sds(spf_1) 0.60 0.49 0.03 1.88 1.00 5289 6287
## sds(sigma_sf0_1) 0.49 0.49 0.01 1.79 1.00 3876 6866
## sds(sigma_spf_1) 0.77 0.61 0.03 2.29 1.00 3279 5300
##
## Group-Level Effects:
## ~dataset (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.47 0.19 0.24 0.95 1.00 4396 6685
## sd(sigma_Intercept) 0.12 0.07 0.02 0.28 1.00 3323 3103
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.13 0.23 -0.32 0.58 1.00 4551 5867
## sigma_Intercept 0.18 0.12 -0.05 0.41 1.00 7990 7681
## sex1 0.16 0.13 -0.08 0.43 1.00 10684 9656
## sigma_sex1 -0.08 0.09 -0.29 0.09 1.00 9147 7879
## sf0_1 -0.28 1.70 -3.76 3.02 1.00 9243 8390
## spf_1 0.71 1.21 -1.85 3.15 1.00 8565 6811
## sigma_sf0_1 -0.28 1.31 -2.61 3.01 1.00 6575 4527
## sigma_spf_1 -1.32 1.73 -5.48 1.38 1.00 5933 6055
## meageage_se -0.01 0.01 -0.02 0.00 1.00 15953 8469
## sigma_meageage_se -0.01 0.00 -0.02 -0.00 1.00 9315 8745
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
| model | Estimate | Est.Error | Q2.5 | Q97.5 | loo_R2 |
|---|---|---|---|---|---|
| h1_after_diagnosis_melsm | 0.1264509 | 0.0160206 | 0.0953042 | 0.1580384 | 0.1098105 |
| h1_after_diagnosis | 0.1266616 | 0.0157463 | 0.0965032 | 0.1587202 | 0.1094989 |
| h1_sex_mod_dataset_varying | 0.1301130 | 0.0154954 | 0.1003210 | 0.1609634 | 0.1085954 |
| h1_simple | 0.1176955 | 0.0153079 | 0.0881099 | 0.1484517 | 0.1045678 |
Participants with lower voice pitch will score lower on neuroticism.
Interpretation: The data are are consistent with linear positive relationship between f0 and neuroticism. The 89% HDI for f0 excludes zero but partly inside the ROPE. 69% of the HDI for pf falls within the ROPE. This evidence is consistent with an association, where people with deeper voices have lower neuroticism (after adjusting for age and gender), but further research is needed to test whether the associations with f0 and pf are negligible in size.
ggplot(vcs, aes(f0, neuro)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).
ggplot(vcs, aes(pf, neuro)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).
ggplot(vcs, aes(age*10, neuro)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
scale_x_continuous("Age") +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing non-finite values (stat_smooth).
## Warning: Removed 796 rows containing missing values (geom_point).
Decisions upon visual diagnosis: Visual diagnostic shows signs of a potential interaction by sex f0 and a potential interaction by sex coupled with nonlinearity for pf. We will fit an interaction for f0, and separate splines by sex for pf. It is necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is included in these analyses.
h1_simple <- brm(neuro ~ f0 + pf + sex + me(age, age_se) + (1 | dataset), data = vcs,
iter = iter, warmup = warmup, chains = chains, cores = chains,
prior = priors,
control = control, save_mevars = TRUE,
file = "models/neuro/h1_simple") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/neuro/h1_simple.rds'
## Automatically saving the model object in 'models/neuro/h1_simple.rds'
## Automatically saving the model object in 'models/neuro/h1_simple.rds'
h1_after_diagnosis <- brm(neuro ~ f0 * sex + s(pf, by = sex) + sex + me(age, age_se) + (1 | dataset), data = vcs,
iter = iter, warmup = warmup, cores = chains, chains = chains,
prior = priors,
control = control, save_mevars = TRUE,
file = "models/neuro/h1_nonlinear") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/neuro/h1_nonlinear.rds'
## Automatically saving the model object in 'models/neuro/h1_nonlinear.rds'
## Automatically saving the model object in 'models/neuro/h1_nonlinear.rds'
conditional_smooths(h1_after_diagnosis)
compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
## elpd_diff se_diff
## h1_after_diagnosis 0.0 0.0
## h1_simple -0.6 2.4
h1_sex_mod_dataset_varying <- brm(neuro ~ f0 + sex + me(age, age_se) +
(1 + f0 + sex + age || dataset), data = vcs,
prior = priors,
iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
control = control, save_mevars = TRUE,
file = "models/neuro/h1_sex_mod_dataset_varying") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/neuro/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/neuro/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/neuro/h1_sex_mod_dataset_varying.rds'
# agree <- brm(
# bf(neuro ~ f0 + sex + me(age, age_se) + (1 | dataset),
# sigma ~ f0 + sex + me(age, age_se) + (1 | dataset)), data = vcs,
# prior = priors,
# iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
# control = control, save_mevars = TRUE,
# file = "models/neuro/agree") %>%
# add_criterion("loo") %>%
# add_criterion("bayes_R2") %>%
# add_criterion("loo_R2")
#
compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying
# ,h1_after_diagnosis_melsm
)
compare_models
## elpd_diff se_diff
## h1_sex_mod_dataset_varying 0.0 0.0
## h1_after_diagnosis -2.1 3.3
## h1_simple -2.7 4.0
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying
# ,h1_after_diagnosis_melsm=h1_after_diagnosis_melsm
) %>%
map(~ { x <- as_tibble(bayes_R2(.))
x$loo_R2 <- loo_R2(.)
x
}) %>%
bind_rows(.id="model") %>%
arrange(desc(loo_R2))
sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
pred.labels = variable_labels)
| Â | Neuroticism | |
|---|---|---|
| Predictors | Estimates | CI (95%) |
| Intercept | 0.04 | -0.28;0.35 |
| Voice pitch (f0) | 0.16 | 0.01;0.30 |
| Formants (f1-f4) | -0.07 | -0.21;0.06 |
| Sex [male] | -0.15 | -0.33;0.02 |
| Age±SE | -0.02 | -0.14;0.09 |
| Random Effects | ||
| σ2 | 0.01 | |
| τ00 | 0.99 | |
| ICC | 0.01 | |
| N dataset | 7 | |
| Observations | 1434 | |
| Marginal R2 / Conditional R2 | 0.088 / 0.088 | |
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
## Possible multicollinearity between b_sex1 and b_f0 (r = 0.77). This might lead to inappropriate results. See 'Details' in '?equivalence_test'.
equiv_test_table <- equiv_test %>%
mutate(Parameter = recode(Parameter, !!!effect_labels)) %>%
mutate(
ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>%
select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
| Parameter | H0 | inside ROPE | 89% HDI |
|---|---|---|---|
| Intercept | Undecided | 51% | [-0.22; 0.30] |
| Voice pitch (f0) | Undecided | 19% | [ 0.03; 0.27] |
| Formants (f1-f4) | Undecided | 69% | [-0.18; 0.04] |
| Sex [male] | Undecided | 25% | [-0.29;-0.01] |
| Age±SE | Undecided | 97% | [-0.11; 0.07] |
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Picking joint bandwidth of 0.00851
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).
Model details
summary(h1_simple)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: neuro ~ f0 + pf + sex + me(age, age_se) + (1 | dataset)
## Data: vcs (Number of observations: 1434)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~dataset (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.14 0.08 0.03 0.34 1.00 1546 1455
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.04 0.16 -0.28 0.35 1.00 2250 3563
## f0 0.16 0.07 0.01 0.30 1.00 4243 5581
## pf -0.07 0.07 -0.21 0.06 1.00 4646 5858
## sex1 -0.15 0.09 -0.33 0.02 1.00 3609 4796
## meageage_se -0.02 0.06 -0.14 0.09 1.00 2512 3192
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.96 0.02 0.92 0.99 1.00 9149 5924
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: neuro ~ f0 * sex + s(pf, by = sex) + sex + me(age, age_se) + (1 | dataset)
## Data: vcs (Number of observations: 1434)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(spfsexfemale_1) 0.63 0.59 0.02 2.16 1.00 4619 3779
## sds(spfsexmale_1) 1.22 0.94 0.07 3.67 1.00 3335 3158
##
## Group-Level Effects:
## ~dataset (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.15 0.08 0.04 0.35 1.00 3066 3700
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.22 0.21 -0.64 0.19 1.00 4524 5308
## f0 0.12 0.08 -0.04 0.27 1.00 7955 6568
## sex1 -0.26 0.13 -0.51 -0.02 1.00 5197 5527
## f0:sex1 -0.16 0.08 -0.31 -0.01 1.00 10869 6230
## spf:sexfemale_1 -0.24 1.50 -3.27 3.03 1.00 4854 4458
## spf:sexmale_1 -0.13 2.07 -4.28 4.26 1.00 6130 5760
## meageage_se -0.01 0.06 -0.12 0.10 1.00 5308 5097
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.95 0.02 0.92 0.99 1.00 10209 5305
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_after_diagnosis_melsm)
# conditional_smooths(h1_after_diagnosis_melsm)
summary(h1_sex_mod_dataset_varying)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: neuro ~ f0 + pf * sex + f0 * sex + me(age, age_se) + (1 + f0 * sex + pf * sex + sex || dataset)
## Data: vcs (Number of observations: 1434)
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~dataset (Number of levels: 7)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.24 0.17 0.02 0.67 1.00 3540 4464
## sd(f0) 0.13 0.12 0.00 0.43 1.00 4402 6207
## sd(sex1) 0.16 0.14 0.01 0.52 1.00 4228 5443
## sd(pf) 0.19 0.17 0.01 0.64 1.00 4928 6626
## sd(f0:sex1) 0.14 0.12 0.01 0.45 1.00 5134 6440
## sd(sex1:pf) 0.31 0.24 0.02 0.89 1.00 3015 4627
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.21 0.22 -0.65 0.21 1.00 9960 8724
## f0 0.13 0.11 -0.08 0.35 1.00 8530 6968
## pf -0.12 0.15 -0.42 0.18 1.00 7696 8006
## sex1 -0.20 0.14 -0.47 0.08 1.00 8816 7551
## pf:sex1 -0.10 0.19 -0.47 0.26 1.00 5858 4908
## f0:sex1 -0.14 0.11 -0.36 0.07 1.00 9590 7033
## meageage_se 0.01 0.06 -0.11 0.12 1.00 11138 8693
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.95 0.02 0.92 0.99 1.00 19230 8270
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# summary(h1_after_diagnosis_melsm)
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
| model | Estimate | Est.Error | Q2.5 | Q97.5 | loo_R2 |
|---|---|---|---|---|---|
| h1_sex_mod_dataset_varying | 0.1065406 | 0.0144428 | 0.0791614 | 0.1352072 | 0.0785213 |
| h1_after_diagnosis | 0.0943398 | 0.0139157 | 0.0682566 | 0.1215960 | 0.0762639 |
| h1_simple | 0.0880062 | 0.0134624 | 0.0622406 | 0.1155778 | 0.0753931 |
Participants with lower voice pitch will report having a more unrestricted sociosexuality.
Interpretation: The data are are consistent with substantial linear negative relationship between f0 and unrestricted sociosexuality. The 89% HDI for f0 falls entirely outside the ROPE, but the HDI for pf falls almost entirely within it. This evidence is consistent with a non-negligible association, where people with deeper voices have a less restricted sociosexuality (after adjusting for age and gender). Further research is needed to verify if the association with pf is truly negligible in size.
ggplot(vcs, aes(f0, soi_full)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing missing values (geom_point).
ggplot(vcs, aes(pf, soi_full)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing missing values (geom_point).
ggplot(vcs, aes(age*10, soi_full)) +
geom_jitter(aes(colour = sex), alpha = 0.3) +
scale_x_continuous("Age") +
geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
geom_smooth(aes(colour = sex, fill = sex), method = "gam",
formula = y ~ s(x)) +
scale_color_viridis_d("Sex") +
scale_fill_viridis_d("Sex")
## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing non-finite values (stat_smooth).
## Warning: Removed 235 rows containing missing values (geom_point).
Decisions upon visual diagnosis: We will fit splines for f0, pf, and age. It seems likely that the nonlinearities in f0 and pf will not both be apparent after adjusting for one another. There is no evidence in the plots that there is an interaction with sex, except for age.
h1_simple <- brm(soi_full ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
iter = iter, warmup = warmup, chains = chains, cores = chains,
prior = priors,
control = control, save_mevars = TRUE,
file = "models/soi_full/h1_simple") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/soi_full/h1_simple.rds'
## Automatically saving the model object in 'models/soi_full/h1_simple.rds'
## Automatically saving the model object in 'models/soi_full/h1_simple.rds'
h1_after_diagnosis <- brm(soi_full ~ s(f0) + s(pf) + sex + s(age, by = sex) + (1 | dataset), data = vcs,
iter = iter, warmup = warmup, cores = chains, chains = chains,
prior = priors,
control = control, save_mevars = TRUE,
file = "models/soi_full/h1_nonlinear") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/soi_full/h1_nonlinear.rds'
## Automatically saving the model object in 'models/soi_full/h1_nonlinear.rds'
## Automatically saving the model object in 'models/soi_full/h1_nonlinear.rds'
conditional_smooths(h1_after_diagnosis)
compare_models <- loo_compare(h1_simple, h1_after_diagnosis)
compare_models
## elpd_diff se_diff
## h1_after_diagnosis 0.0 0.0
## h1_simple -3.5 3.6
h1_sex_mod_dataset_varying <- brm(soi_full ~ f0 + sex + s(age, by = sex) +
(1 + f0 + sex + age || dataset), data = vcs,
prior = priors,
iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
control = control, save_mevars = TRUE,
file = "models/soi_full/h1_sex_mod_dataset_varying") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/soi_full/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/soi_full/h1_sex_mod_dataset_varying.rds'
## Automatically saving the model object in 'models/soi_full/h1_sex_mod_dataset_varying.rds'
h1_after_diagnosis_melsm <- brm(
bf(soi_full ~ f0 + sex + s(age, by = sex) + (1 | dataset),
sigma ~ f0 + sex + s(age, by = sex) + (1 | dataset)), data = vcs,
prior = priors,
iter = iter + 2000, warmup = warmup + 1000, cores = chains, chains = chains,
control = control, save_mevars = TRUE,
file = "models/soi_full/h1_after_diagnosis_melsm") %>%
add_criterion("loo") %>%
add_criterion("bayes_R2") %>%
add_criterion("loo_R2")
## Automatically saving the model object in 'models/soi_full/h1_after_diagnosis_melsm.rds'
## Automatically saving the model object in 'models/soi_full/h1_after_diagnosis_melsm.rds'
## Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
## Automatically saving the model object in 'models/soi_full/h1_after_diagnosis_melsm.rds'
compare_models <- loo_compare(h1_simple, h1_after_diagnosis, h1_sex_mod_dataset_varying, h1_after_diagnosis_melsm)
compare_models
## elpd_diff se_diff
## h1_after_diagnosis_melsm 0.0 0.0
## h1_after_diagnosis -8.2 5.5
## h1_sex_mod_dataset_varying -11.4 5.4
## h1_simple -11.7 6.0
r2s <- list(h1_simple=h1_simple, h1_after_diagnosis=h1_after_diagnosis, h1_sex_mod_dataset_varying=h1_sex_mod_dataset_varying, h1_after_diagnosis_melsm=h1_after_diagnosis_melsm) %>%
map(~ { x <- as_tibble(bayes_R2(.))
x$loo_R2 <- loo_R2(.)
x
}) %>%
bind_rows(.id="model") %>%
arrange(desc(loo_R2))
sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
pred.labels = variable_labels)
| Â |
Unrestricted sociosexuality |
|
|---|---|---|
| Predictors | Estimates | CI (95%) |
| Intercept | -0.20 | -0.60;0.19 |
| Voice pitch (f0) | -0.29 | -0.41;-0.17 |
| Formants (f1-f4) | -0.01 | -0.14;0.12 |
| Sex [male] | -0.03 | -0.17;0.12 |
| Age | 0.09 | -0.01;0.18 |
| Random Effects | ||
| σ2 | 0.08 | |
| τ00 | 0.93 | |
| ICC | 0.07 | |
| N dataset | 9 | |
| Observations | 1995 | |
| Marginal R2 / Conditional R2 | 0.160 / 0.160 | |
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>%
mutate(Parameter = recode(Parameter, !!!effect_labels)) %>%
mutate(
ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>%
select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")
| Parameter | H0 | inside ROPE | 89% HDI |
|---|---|---|---|
| Intercept | Undecided | 26% | [-0.51; 0.13] |
| Voice pitch (f0) | Rejected | 0% | [-0.39;-0.19] |
| Formants (f1-f4) | Undecided | 96% | [-0.12; 0.09] |
| Sex [male] | Undecided | 89% | [-0.14; 0.09] |
| Age | Undecided | 60% | [ 0.02; 0.16] |
plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Picking joint bandwidth of 0.0073
## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).
Model details
summary(h1_simple)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: soi_full ~ f0 + pf + sex + age + (1 | dataset)
## Data: vcs (Number of observations: 1995)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~dataset (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.46 0.15 0.27 0.85 1.00 1910 3239
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.20 0.20 -0.60 0.19 1.00 2156 3365
## f0 -0.29 0.06 -0.41 -0.17 1.00 4808 4766
## pf -0.01 0.06 -0.14 0.12 1.00 6259 5195
## sex1 -0.03 0.07 -0.17 0.12 1.00 4258 4948
## age 0.09 0.05 -0.01 0.18 1.00 7240 5169
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.92 0.01 0.89 0.95 1.00 7124 4467
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: soi_full ~ s(f0) + s(pf) + sex + s(age, by = sex) + (1 | dataset)
## Data: vcs (Number of observations: 1995)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sds(sf0_1) 1.62 1.08 0.08 4.11 1.00 1484
## sds(spf_1) 0.67 0.58 0.02 2.15 1.00 3178
## sds(sagesexfemale_1) 1.06 0.68 0.19 2.80 1.00 2965
## sds(sagesexmale_1) 0.49 0.46 0.02 1.72 1.00 3657
## Tail_ESS
## sds(sf0_1) 2166
## sds(spf_1) 4334
## sds(sagesexfemale_1) 2570
## sds(sagesexmale_1) 4202
##
## Group-Level Effects:
## ~dataset (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.44 0.14 0.25 0.80 1.00 2878 4120
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.02 0.16 -0.31 0.33 1.00 1905 3744
## sex1 0.02 0.12 -0.20 0.26 1.00 6673 6118
## sf0_1 -1.68 2.06 -5.25 3.00 1.00 4230 5063
## spf_1 0.02 1.26 -2.52 2.75 1.00 6463 5365
## sage:sexfemale_1 0.88 1.70 -2.67 4.38 1.00 5886 5216
## sage:sexmale_1 0.40 1.05 -1.82 2.67 1.00 5799 4865
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.92 0.01 0.89 0.95 1.00 12164 5749
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(h1_simple)
conditional_smooths(h1_after_diagnosis)
summary(h1_sex_mod_dataset_varying)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: soi_full ~ f0 + pf * sex + s(age, by = sex) + (1 + pf + f0 + sex || dataset)
## Data: vcs (Number of observations: 1995)
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
## total post-warmup samples = 12000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sds(sagesexfemale_1) 1.01 0.66 0.17 2.68 1.00 5601
## sds(sagesexmale_1) 0.62 0.55 0.03 2.04 1.00 4981
## Tail_ESS
## sds(sagesexfemale_1) 4654
## sds(sagesexmale_1) 6139
##
## Group-Level Effects:
## ~dataset (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.45 0.16 0.24 0.84 1.00 5523 7210
## sd(pf) 0.13 0.11 0.01 0.39 1.00 4931 6390
## sd(f0) 0.10 0.08 0.00 0.30 1.00 4961 6208
## sd(sex1) 0.13 0.13 0.00 0.47 1.00 4191 5638
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.03 0.17 -0.39 0.31 1.00 4576 6173
## f0 -0.29 0.08 -0.44 -0.14 1.00 10385 7627
## pf -0.02 0.10 -0.22 0.18 1.00 9863 8486
## sex1 -0.02 0.12 -0.26 0.21 1.00 8065 7178
## pf:sex1 -0.03 0.09 -0.21 0.14 1.00 11580 8726
## sage:sexfemale_1 0.70 1.67 -2.84 4.02 1.00 8781 8376
## sage:sexmale_1 0.62 1.22 -1.83 3.28 1.00 7634 6686
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.92 0.01 0.89 0.95 1.00 18769 8603
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(h1_after_diagnosis_melsm)
## Warning: There were 10 divergent transitions after warmup. Increasing
## adapt_delta above 0.99 may help. See http://mc-stan.org/misc/
## warnings.html#divergent-transitions-after-warmup
## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: soi_full ~ f0 + pf + sex + s(age, by = sex) + (1 | dataset)
## sigma ~ f0 + pf + sex + s(age, by = sex) + (1 | dataset)
## Data: vcs (Number of observations: 1995)
## Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
## total post-warmup samples = 12000
##
## Smooth Terms:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sds(sagesexfemale_1) 1.11 0.67 0.26 2.79 1.00 6683
## sds(sagesexmale_1) 0.44 0.42 0.01 1.58 1.00 6492
## sds(sigma_sagesexfemale_1) 0.63 0.61 0.02 2.25 1.00 3904
## sds(sigma_sagesexmale_1) 0.40 0.43 0.01 1.59 1.00 5024
## Tail_ESS
## sds(sagesexfemale_1) 6966
## sds(sagesexmale_1) 7184
## sds(sigma_sagesexfemale_1) 7777
## sds(sigma_sagesexmale_1) 6947
##
## Group-Level Effects:
## ~dataset (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.45 0.15 0.26 0.82 1.00 4973 7393
## sd(sigma_Intercept) 0.17 0.07 0.08 0.33 1.00 3835 6626
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 0.02 0.16 -0.29 0.35 1.00 3550
## sigma_Intercept -0.12 0.06 -0.25 0.01 1.00 4597
## f0 -0.29 0.06 -0.41 -0.17 1.00 16252
## pf 0.01 0.07 -0.12 0.14 1.00 18032
## sex1 -0.00 0.07 -0.15 0.14 1.00 12734
## sigma_f0 -0.00 0.05 -0.09 0.09 1.00 16967
## sigma_pf 0.03 0.05 -0.08 0.13 1.00 17853
## sigma_sex1 0.05 0.06 -0.06 0.16 1.00 12697
## sage:sexfemale_1 0.80 1.74 -2.91 4.29 1.00 10802
## sage:sexmale_1 0.25 0.98 -1.86 2.29 1.00 9187
## sigma_sage:sexfemale_1 -0.46 1.75 -4.77 2.73 1.00 7930
## sigma_sage:sexmale_1 -0.11 1.03 -2.06 2.32 1.00 7920
## Tail_ESS
## Intercept 4945
## sigma_Intercept 5941
## f0 8933
## pf 9395
## sex1 8948
## sigma_f0 9549
## sigma_pf 9052
## sigma_sex1 8763
## sage:sexfemale_1 8699
## sage:sexmale_1 7383
## sigma_sage:sexfemale_1 5720
## sigma_sage:sexmale_1 5717
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
kable(r2s, caption = "Explained variance, regular and leave-one-out.")
| model | Estimate | Est.Error | Q2.5 | Q97.5 | loo_R2 |
|---|---|---|---|---|---|
| h1_after_diagnosis | 0.1688764 | 0.0139599 | 0.1416838 | 0.1963915 | 0.1537734 |
| h1_after_diagnosis_melsm | 0.1679133 | 0.0121578 | 0.1442473 | 0.1917065 | 0.1528378 |
| h1_sex_mod_dataset_varying | 0.1688764 | 0.0140204 | 0.1418163 | 0.1967623 | 0.1512160 |
| h1_simple | 0.1598969 | 0.0137611 | 0.1322719 | 0.1866512 | 0.1509295 |